Obtaining smooth solutions to large, linear, inverse problems

It is not uncommon now for geophysical inverse problems to be parameterized by 10 to 105 unknowns associated with upwards of 106 to 107 data constraints. The matrix problem defining the linearization of such a system (e.g., ~m = b) is usually solved with a least-squares criterion (m = (~t~ -) ~ tb). The size of the matrix, however, discourages the direct solution of the system and researchers often tum to iterative techniques such as the method of conjugate gradients to obtain an estimate of the least-squares solution. These iterative methods take advantage of the sparseness of ~, which often has as few as 2-3 percent of its elements nonzero, and do not require the calculation (or storage) of the matrix ~t~. Although there are usually many more data constraints than unknowns, these problems are, in general, underdetermined and therefore require some sort of regularization to obtain a solution. When the regularization is simple damping, the conjugate gradients method tends to converge in relativelyfew iterations. However, when derivative-type regularization is applied (first derivative constraints to obtain the flattest model that fits the data; second derivative to obtain the smoothest), the convergence of parts ofthe solution may be drastically inhibited. In a series of I-D examples and a synthetic 2-D crosshole tomography example, we demonstrate this problem and also suggest a method of accelerating the convergence through the preconditioning of the conjugate gradient search directions. We derive a I-D preconditioning operator for the case of first derivative regularization using a WKBJ approximation. We have found that preconditioning can reduce the number of iterations necessary to obtain satisfactory convergence by up to an order of magnitude. The conclusions we present are also relevant to Bayesian inversion, where a smoothness constraint is imposed through an a priori covariance of the model.